/*******************************************************************************

				HDFUTIL
  Title:        I/O functions for writing to or reading from HDF fies..    
  Author:       Hanna Szoke
  Date:         12/04/05
  Function:     This package provides all interfaces between Speden and H5F.  

*******************************************************************************/

#include "speden.h"	/* includes util.h, cellparams.h, symmetry.h, dims.h */
#include "sp.h"
#include <hdf5.h>

hid_t   file_id;           
hid_t   group_id;
hid_t   dataspace_id;
hid_t   dataset_id;

herr_t  status;
herr_t  H5Eset_auto(H5E_auto_t, void *) ;

hsize_t dims[3];
static     char	ufilename[MAXSTRING] ;

#ifdef DOUBLE
	#define MEM_DATA_TYPE  H5T_NATIVE_DOUBLE
#else
	#define MEM_DATA_TYPE  H5T_NATIVE_FLOAT
#endif

        /************************************************************
        Read a file of electrons per voxel in HDF format.
	The file is expected to exist; otherwise, code fails fatally.
        ************************************************************/
void    read_HDFnump(char *filename, real *nump)
{
     H5E_auto_t	H5Eprint = 0 ;

     status = H5Eset_auto(NULL, NULL) ;  /* toggle OFF */
     file_id = H5Fopen(filename, H5F_ACC_RDWR, H5P_DEFAULT);
     if (file_id < 0) {
          sprintf(message, "Couldn't open %s\n", filename) ;
          Fatal(message) ;
     }

     status = H5Eset_auto(H5Eprint, stderr) ;  /* toggle ON */

     dataset_id = H5Dopen(file_id, "/nump");
     status = H5Dread(dataset_id, MEM_DATA_TYPE, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, nump);
     status = H5Dclose(dataset_id);

     return ;
}

        /************************************************************
        Get a unique name for writing an HDF file.
        ************************************************************/
void get_unique_hdfname(char *inname)
{
     int	j ;

     H5E_auto_t	H5Eprint = 0 ;

     sprintf(ufilename, "%s.hdf", inname) ;

     status = H5Eset_auto(NULL, NULL) ;  /* toggle OFF */
     file_id = H5Fopen(ufilename, H5F_ACC_RDWR, H5P_DEFAULT);
     status = H5Eset_auto(H5Eprint, stderr) ;  /* toggle ON */

     if (file_id >= 0){
    
     /* Basic name is in use; find first free name of form infile[j]  */

         H5Fclose(file_id);

         for (j = 1; j < 10; j++) {

             sprintf(ufilename, "%s%1d.hdf", inname, j) ;

             status = H5Eset_auto(NULL, NULL) ;  /* toggle OFF */
             file_id = H5Fopen(ufilename, H5F_ACC_RDWR, H5P_DEFAULT);
             status = H5Eset_auto(H5Eprint, stderr) ;  /* toggle ON */
             if (file_id < 0)
                  break ;
             else
                  H5Fclose(file_id) ;
         }

         if (j == 9) 
              Fatal("\nThere are no more free names for the hdf file!\n") ;
     }
     sprintf(message, "\nWriting hdf data to file %s", ufilename) ;
     printTwice(message) ;

     return ;
}

        /************************************************************
        Write a file of electrons per voxel in HDF format.
        ************************************************************/

void    write_HDFnump(char *filename, real *nump)
{
     void	get_unique_hdfname(char *) ;

     get_unique_hdfname(filename) ;

     file_id = H5Fcreate(ufilename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
     group_id = H5Gcreate(file_id, "/nump", 0);

     dims[0] = Nx;
     dims[1] = Ny;
     dims[2] = Nz;

     dataspace_id = H5Screate_simple( 3, dims, NULL );
     
     dataset_id = H5Dcreate(file_id, "/nump/data", H5T_NATIVE_FLOAT,
                  dataspace_id, H5P_DEFAULT);

     status = H5Dwrite(dataset_id, MEM_DATA_TYPE, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, nump);

     status = H5Dclose(dataset_id);
     status = H5Sclose(dataspace_id);

     status = H5Gclose(group_id);
     status = H5Fclose(file_id);

     return ;
}

        /************************************************************
        Read a file of electrons densities in HDF format.
	The file is expected to exist; otherwise, code fails fatally.
        ************************************************************/

void    read_HDFmap(char *filename, real *din, real ar_min, real ar_max)
{
     H5E_auto_t	H5Eprint = 0 ;
     float 	values[2] ;

     status = H5Eset_auto(NULL, NULL) ;  /* toggle OFF */
     file_id = H5Fopen(filename, H5F_ACC_RDWR, H5P_DEFAULT);
     if (file_id < 0) {
          sprintf(message, "Couldn't open %s\n", filename) ;
          Fatal(message) ;
     }

     status = H5Eset_auto(H5Eprint, stderr) ;  /* toggle ON */

     dataset_id = H5Dopen(file_id, "/map/data_range") ;
     status = H5Dread(dataset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, values);
     status = H5Dclose(dataset_id);

     ar_min = values[0] ;
     ar_max = values[1] ;
     
     dataset_id = H5Dopen(file_id, "/map/data") ;
     status = H5Dread(dataset_id, MEM_DATA_TYPE, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, din);
     status = H5Dclose(dataset_id);

     status = H5Fclose(file_id);

     return ;
}
	/************************************************************
	Write a file of output electron densities in HDF format.
	If file existed, it is destroyed (quietly).
	************************************************************/

void	write_HDFmap(char *filename, real *dout, real ar_min, real ar_max)	
{
     float 	values[2] ;
     void	get_unique_hdfname(char *) ;

     get_unique_hdfname(filename) ;

     file_id = H5Fcreate(ufilename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT) ;
     group_id = H5Gcreate(file_id, "/map", 0);

     dims[0] = 2 ;
     dataspace_id = H5Screate_simple( 1, dims, NULL ); 
     dataset_id = H5Dcreate(file_id, "/map/data_range", H5T_NATIVE_FLOAT, 
                  dataspace_id, H5P_DEFAULT);
     values[0] = ar_min ;
     values[1] = ar_max ;
     status = H5Dwrite(dataset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, 
              H5P_DEFAULT, values);
     status = H5Dclose(dataset_id);
     status = H5Sclose(dataspace_id);

     dims[0] = Nx;
     dims[1] = Ny;
     dims[2] = Nz;
     dataspace_id = H5Screate_simple( 3, dims, NULL ); 

     dataset_id = H5Dcreate(file_id, "/map/data", H5T_NATIVE_FLOAT, 
                  dataspace_id, H5P_DEFAULT);

     status = H5Dwrite(dataset_id, MEM_DATA_TYPE, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, dout);

        /* End access to the dataset and release resources used by it. */
     status = H5Dclose(dataset_id);
     status = H5Sclose(dataspace_id);

        /* Close the group and file */
     status = H5Gclose(group_id);

     status = H5Fclose( file_id );

     return ;
}

        /************************************************************
        Read data length, N, from a file in HDF format.  The file 
	and data are expected to exist; otherwise, code fails fatally.
        ************************************************************/

int	HDFdatalength(char *filename, char *dset)
{
     H5E_auto_t H5Eprint = 0 ;
     int	N ;

     status = H5Eset_auto(NULL, NULL) ;  /* toggle OFF */
     file_id = H5Fopen(filename, H5F_ACC_RDWR, H5P_DEFAULT);
     if (file_id < 0) {
          sprintf(message, "Couldn't open %s\n", filename) ;
          Fatal(message) ;
     }

     dataset_id = H5Dopen(file_id, dset) ;
     if (dataset_id < 0) {
          sprintf(message, "Couldn't find %s in %s\n", dset, filename) ;
          Fatal(message) ;
     }

     status = H5Eset_auto(H5Eprint, stderr) ;  /* toggle ON */
     status = H5Dread(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, &N);
     status = H5Dclose(dataset_id);
     status = H5Fclose(file_id);

     return (N) ;

}
        /************************************************************
        Read a file of fcalc information in HDF format.
        The file is expected to exist; otherwise, code fails fatally.
        ************************************************************/

int    read_HDFfcalc(char *filename) 
{
     int	compress_fc() ;
     H5E_auto_t H5Eprint = 0 ;

     status = H5Eset_auto(NULL, NULL) ;  /* toggle OFF */
     file_id = H5Fopen(filename, H5F_ACC_RDWR, H5P_DEFAULT);
     if (file_id < 0) {
          sprintf(message, "Couldn't open %s\n", filename) ;
          Fatal(message) ;
     }
     status = H5Eset_auto(H5Eprint, stderr) ;  /* toggle ON */


     dataset_id = H5Dopen(file_id, "/fcalc/h") ;
     status = H5Dread(dataset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Fc_Hlist);
     status = H5Dclose(dataset_id);

     dataset_id = H5Dopen(file_id, "/fcalc/k") ;
     status = H5Dread(dataset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Fc_Klist);
     status = H5Dclose(dataset_id);

     dataset_id = H5Dopen(file_id, "/fcalc/l") ;
     status = H5Dread(dataset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Fc_Llist);
     status = H5Dclose(dataset_id);

     dataset_id = H5Dopen(file_id, "/fcalc/re") ;
     status = H5Dread(dataset_id, MEM_DATA_TYPE, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Fc_relist);
     status = H5Dclose(dataset_id);

     dataset_id = H5Dopen(file_id, "/fcalc/im") ;
     status = H5Dread(dataset_id, MEM_DATA_TYPE, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Fc_imlist);
     status = H5Dclose(dataset_id);

     status = H5Fclose(file_id);

     Nfc = compress_fc() ;

     return (Nfc);
}

        /************************************************************
	Compress entries based on current problem's resolution;
	This is the equivalent of the old set_off_grid_fc().
        ************************************************************/
int	compress_fc()
{
     int	nin, nout ;
     int	out ;
     float	h, k, l ;
     real	re0, im0 ;

     for (nin = nout = 0; nin < Nfob; nin++) {
         out = FALSE ;
         h = *(Fc_Hlist+nin) ;
         k = *(Fc_Klist+nin) ;
         l = *(Fc_Llist+nin) ;

         if ((h==0)&&(k==0)&&(l==0)) {
             re0 = *(Fc_relist+nin) ;
             im0 = *(Fc_imlist+nin) ;
             amp0 = sqrt(re0*re0 + im0*im0) ;
             if (nin > 0) {
                  sprintf (message,
                      "F(000) is the %d-th entry in your input.", nin) ;
                  Warning(message) ;
             }
         }
         if (h < 0) {
             h = -h ;
             k = -k ;
             l = -l ;
         }
         if (h > Nh-1)
             out = TRUE ;
         if ((k < -Nk/2) || (k >= Nk/2))
             out = TRUE ;
         if ((l < -Nl/2) || (l >= Nl/2))
             out = TRUE ;

         if (!out) {
              *(Fc_Hlist + nout) = h ;
              *(Fc_Klist + nout) = k ;
              *(Fc_Llist + nout) = l ;
              *(Fc_relist + nout) = *(Fc_relist + nin) ;
              *(Fc_imlist + nout) = *(Fc_imlist + nin) ;
              nout++ ;
         }
     }
     return(nout) ;
}

        /************************************************************
        Write a file of output fcalcs in HDF format.
        ************************************************************/

void    write_HDFfcalc(char *filename)
{
     void	get_unique_hdfname(char *) ;
     float	values[1] ;	

     get_unique_hdfname(filename) ;
     file_id = H5Fcreate(ufilename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
     group_id = H5Gcreate(file_id, "/fcalc", 0);

     dims[0] = 1 ;
     dataspace_id = H5Screate_simple( 1, dims, NULL );
     dataset_id = H5Dcreate(file_id, "/fcalc/Nfc", H5T_NATIVE_FLOAT,
                  dataspace_id, H5P_DEFAULT);
     values[0] = Nfc ;
     status = H5Dwrite(dataset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, values);
     status = H5Dclose(dataset_id);
     status = H5Sclose(dataspace_id);

     dims[0] = Nfc;
     dataspace_id = H5Screate_simple( 1, dims, NULL );
     dataset_id = H5Dcreate(file_id, "/fcalc/h", H5T_NATIVE_FLOAT,
                  dataspace_id, H5P_DEFAULT);
     status = H5Dwrite(dataset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Fc_Hlist);
     status = H5Dclose(dataset_id);
     status = H5Sclose(dataspace_id);

     dataspace_id = H5Screate_simple( 1, dims, NULL );
     dataset_id = H5Dcreate(file_id, "/fcalc/k", H5T_NATIVE_FLOAT,
                  dataspace_id, H5P_DEFAULT);
     status = H5Dwrite(dataset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Fc_Klist);
     status = H5Dclose(dataset_id);
     status = H5Sclose(dataspace_id);

     dataspace_id = H5Screate_simple( 1, dims, NULL );
     dataset_id = H5Dcreate(file_id, "/fcalc/l", H5T_NATIVE_FLOAT,
                  dataspace_id, H5P_DEFAULT);
     status = H5Dwrite(dataset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Fc_Llist);
     status = H5Dclose(dataset_id);
     status = H5Sclose(dataspace_id);

     dataspace_id = H5Screate_simple( 1, dims, NULL );
     dataset_id = H5Dcreate(file_id, "/fcalc/re", H5T_NATIVE_FLOAT,
                  dataspace_id, H5P_DEFAULT);
     status = H5Dwrite(dataset_id, MEM_DATA_TYPE, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Fc_relist);
     status = H5Dclose(dataset_id);
     status = H5Sclose(dataspace_id);

     dataspace_id = H5Screate_simple( 1, dims, NULL );
     dataset_id = H5Dcreate(file_id, "/fcalc/im", H5T_NATIVE_FLOAT,
                  dataspace_id, H5P_DEFAULT);
     status = H5Dwrite(dataset_id, MEM_DATA_TYPE, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Fc_imlist);
     status = H5Dclose(dataset_id);
     status = H5Sclose(dataspace_id);

        /* Close the group and file */
     status = H5Gclose(group_id);
     status = H5Fclose( file_id );

     return ;
}

        /************************************************************
        Read a file of fobs information in HDF format.
        The file is expected to exist; otherwise, code fails fatally.
        ************************************************************/

int    read_HDFfobs(char *filename) 
{
     int	compress_fo() ;
     H5E_auto_t H5Eprint = 0 ;

     status = H5Eset_auto(NULL, NULL) ;  /* toggle OFF */
     file_id = H5Fopen(filename, H5F_ACC_RDWR, H5P_DEFAULT);
     if (file_id < 0) {
          sprintf(message, "Couldn't open %s\n", filename) ;
          Fatal(message) ;
     }
     status = H5Eset_auto(H5Eprint, stderr) ;  /* toggle ON */

     dataset_id = H5Dopen(file_id, "/fobs/Nfob") ;
     status = H5Dread(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, &Nfob);
     status = H5Dclose(dataset_id);

     dataset_id = H5Dopen(file_id, "/fobs/h") ;
     status = H5Dread(dataset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Hlist);
     status = H5Dclose(dataset_id);

     dataset_id = H5Dopen(file_id, "/fobs/k") ;
     status = H5Dread(dataset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Klist);
     status = H5Dclose(dataset_id);

     dataset_id = H5Dopen(file_id, "/fobs/l") ;
     status = H5Dread(dataset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Llist);
     status = H5Dclose(dataset_id);

     dataset_id = H5Dopen(file_id, "/fobs/amp") ;
     status = H5Dread(dataset_id, MEM_DATA_TYPE, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Folist);
     status = H5Dclose(dataset_id);

     dataset_id = H5Dopen(file_id, "/fobs/sig") ;
     status = H5Dread(dataset_id, MEM_DATA_TYPE, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Siglist);
     status = H5Dclose(dataset_id);

     status = H5Fclose(file_id);

     /*	Compress entries based on current problem's resolution	*/

     Nfob = compress_fo() ;

     return (Nfob) ;
}

        /************************************************************
	Compress entries based on current problem's resolution;
	This is the equivalent of the old "set_fo()".
        ************************************************************/
int	compress_fo()
{
     int	nin, nout ;
     real	fo_nusq(int) ;
     int	out ;
     float	h, k, l ;

     for (nin = nout = 0; nin < Nfob; nin++) {
         out = FALSE ;
         h = *(Hlist+nin) ;
         k = *(Klist+nin) ;
         l = *(Llist+nin) ;

          if ((h==0)&&(k==0)&&(l==0)) {
              if (nin == 0) {
                 F000 = *(Folist) ;
                 S000 = *(Siglist) ;
              }
              else {
                  sprintf (message,
                 "F(000) is the %d-th entry in your input.", nin) ;
                  Warning(message) ;
		  /*                  Fatal("Quitting") ;*/
              }
         }
/***
         if (h < 0) {
             h = -h ;
             k = -k ;
             l = -l ;
         }
***/
         if (h > Nh-1)
             out = TRUE ;
         if ((k < -Nk/2) || (k >= Nk/2))
             out = TRUE ;
         if ((l < -Nl/2) || (l >= Nl/2))
             out = TRUE ;
         if (fo_nusq(nin) >= limit_rat) 
             out = TRUE ;

         if (!out) {
              *(Hlist + nout) = h ;
              *(Klist + nout) = k ;
              *(Llist + nout) = l ;
              *(Folist  + nout) = *(Folist  + nin) ;
              *(Siglist + nout) = *(Siglist + nin) ;
              nout++ ;
         }
     }
     return(nout) ;
}

        /************************************************************
        Write a file of output fobs in HDF format.
        ************************************************************/

void    write_HDFfobs(char *filename)
{
     void	get_unique_hdfname(char *) ;
     float	values[1] ;	

     get_unique_hdfname(filename) ;
     file_id = H5Fcreate(ufilename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
     group_id = H5Gcreate(file_id, "/fobs", 0);

     dims[0] = 1 ;
     dataspace_id = H5Screate_simple( 1, dims, NULL );
     dataset_id = H5Dcreate(file_id, "/fobs/Nfob", H5T_NATIVE_INT,
                  dataspace_id, H5P_DEFAULT);
     values[0] = Nfob ;
     status = H5Dwrite(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, values);
     status = H5Dclose(dataset_id);
     status = H5Sclose(dataspace_id);

     dims[0] = Nfob;
     dataspace_id = H5Screate_simple( 1, dims, NULL );
     dataset_id = H5Dcreate(file_id, "/fobs/h", H5T_NATIVE_FLOAT,
                  dataspace_id, H5P_DEFAULT);
     status = H5Dwrite(dataset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Fc_Hlist);
     status = H5Dclose(dataset_id);
     status = H5Sclose(dataspace_id);

     dataspace_id = H5Screate_simple( 1, dims, NULL );
     dataset_id = H5Dcreate(file_id, "/fobs/k", H5T_NATIVE_FLOAT,
                  dataspace_id, H5P_DEFAULT);
     status = H5Dwrite(dataset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Fc_Klist);
     status = H5Dclose(dataset_id);
     status = H5Sclose(dataspace_id);

     dataspace_id = H5Screate_simple( 1, dims, NULL );
     dataset_id = H5Dcreate(file_id, "/fobs/l", H5T_NATIVE_FLOAT,
                  dataspace_id, H5P_DEFAULT);
     status = H5Dwrite(dataset_id, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Fc_Llist);
     status = H5Dclose(dataset_id);
     status = H5Sclose(dataspace_id);

     dataspace_id = H5Screate_simple( 1, dims, NULL );
     dataset_id = H5Dcreate(file_id, "/fobs/amp", H5T_NATIVE_FLOAT,
                  dataspace_id, H5P_DEFAULT);
     status = H5Dwrite(dataset_id, MEM_DATA_TYPE, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Fc_relist);
     status = H5Dclose(dataset_id);
     status = H5Sclose(dataspace_id);

     dataspace_id = H5Screate_simple( 1, dims, NULL );
     dataset_id = H5Dcreate(file_id, "/fobs/sig", H5T_NATIVE_FLOAT,
                  dataspace_id, H5P_DEFAULT);
     status = H5Dwrite(dataset_id, MEM_DATA_TYPE, H5S_ALL, H5S_ALL,
              H5P_DEFAULT, Fc_imlist);
     status = H5Dclose(dataset_id);
     status = H5Sclose(dataspace_id);

        /* Close the group and file */
     status = H5Gclose(group_id);
     status = H5Fclose( file_id );

     return ;
}
